home *** CD-ROM | disk | FTP | other *** search
/ EnigmA Amiga Run 1997 July / EnigmA AMIGA RUN 20 (1997)(G.R. Edizioni)(IT)[!][issue 1997-07 & 08][EAR-CD IV].iso / earcd / dev / c / amcsrc2.lha / AMCSources2 / FFT / IFFT.c < prev    next >
C/C++ Source or Header  |  1996-07-17  |  9KB  |  389 lines
  1. /*    name...
  2.         ifft
  3.  
  4.     bugs...
  5.         Should do transforms on several functions (-b option).
  6.  
  7.     history...
  8.         -- 1.21 --
  9.         1 Jun 87    -a option suppresses abscissas on output.
  10.         -- 1.20 --
  11.         11 May 87    Using "float" rather than "double" variables
  12.         8 May 87    Reading into doubles and performing range check
  13.                     before converting to floats.
  14.         -- 1.1 --
  15.         6 Apr 87    Default abscissa step is 1 (rather than 0!)
  16.                     Some printing removed.
  17.         -- 1.0 --
  18.         20 Mar 87    default output to STDOUT
  19.         8 Feb 87    split off from fft program.
  20.  
  21.     performance...
  22.         128 point transform in 14.43 sec with 7.5 MHz V-20 and no 8087.
  23.         4096 point transform in 1:53.15 with 6 MHz 80286 and 4 MHz 80287.
  24. */
  25. #include <stdio.h>
  26. #include <math.h>
  27.  
  28. #define VERSION "1.21"
  29.  
  30.  
  31. #define BIGFLOAT 6.e38    /* no floats larger than this  */
  32. #define ENTRIES 4098
  33. #define MAXLABELS 50
  34. #define BUFSIZE 200
  35. #define FLOAT float    /* change to "double" to restore higher precision */
  36.  
  37. FILE *ifile=stdin, *ofile=stdout;
  38.  
  39. char buffer[128], iname[35], oname[35];
  40. char buf[BUFSIZE];
  41. char *label_array[MAXLABELS];
  42. char **p_text=label_array;
  43.  
  44.  
  45. int m,            /* number of input values to be retained */
  46. n,                /* number of data points to be Fourier transformed */
  47. *nnn,            /* points to n */
  48. breaking,        /* nonzero if finding separate transforms for several functions */
  49. labels,            /* number of labels stored */
  50. automatic_abscissas, /* nonzero if abscissas to be calculated */
  51. abscissa_arguments,
  52. keeping=0,        /* if nonzero, the number of input values to be kept */
  53. f_arguments=0,
  54. x_arguments=0,
  55. index_array[MAXLABELS],    /* indexes into x and y */
  56. *p_data=index_array;
  57.  
  58. FLOAT *x,        /* frequency values */
  59. *y;                /* magnitudes (on input) or voltages (on output) */
  60. double abscissa,
  61. origin,            /* abscissa value to be treated as time zero */
  62. abscissa_step=1.,
  63. fmin, fmax,
  64. xmin, xmax;        /* first & last data points to be used */
  65.  
  66. main(int argc,char **argv)
  67. {    int i;
  68.  
  69.     {double junk;
  70.     sscanf("3.4","%lf",&junk);        /* workaround for DESMET sscanf bug */
  71.     }
  72.  
  73.     argc = args(argc, argv);        /* fetch switches */
  74.  
  75.     read_data(argc, argv);
  76.     pad();
  77.  
  78. /*    fprintf(stderr, "\n\n   computing inverse fft\n\n");  */
  79.     ifft(&n, x, y);
  80.  
  81.     i=1;
  82.     while(1)
  83.         {if(!automatic_abscissas)  fprintf(ofile, "%15.6g", x[i] + origin);
  84.         fprintf(ofile, "%15.6g", y[i]);
  85.         if(++i>n) break;
  86.         fprintf(ofile, "\n");
  87.         }
  88.     fprintf(ofile, " \"\"\n");
  89.     if(ofile!=stdout)
  90.         {fclose(ofile);
  91.         }
  92. }
  93.  
  94. pad()
  95. {    int m;
  96.     double freq, df;
  97.     for (m=2; m<ENTRIES; m *= 2)
  98.         if (n<=m+1) break;
  99.     if(n==m+1) return;
  100. /*    printf("transforming %d points after padding\n", m+1);  */
  101.     freq=x[n];
  102.     df=freq-x[n-1];
  103.     while(n<=m)
  104.         {n++;
  105.         freq += df;
  106.         x[n]=freq;
  107.         y[2*n-1]=y[2*n]=0.;
  108.         }
  109. }
  110.  
  111. listt(y, i, j) FLOAT y[]; int i, j;    /* display y[i] through y[j]  */
  112. {    while (i<=j) {printf("%4d %15.9f \n", i, y[i]); ++i;}
  113. }
  114.  
  115.  
  116. read_data(argc, argv) int argc; char **argv;
  117. {    int i;
  118. //        int j;
  119.         int length, nums;
  120.     double xx, yyr, yyi;
  121. //        double d, *pd, sum;
  122.     char *s, *t;
  123. //        char *stop;
  124.     char *malloc();
  125.     char *strchr();
  126.  
  127.     x=(FLOAT *)malloc(ENTRIES*sizeof(FLOAT));
  128.     y=(FLOAT *)malloc(2*ENTRIES*sizeof(FLOAT));
  129.     if(x==0 || y==0) {fprintf(stderr, "can\'t allocate buffer"); exit(1);}
  130.     argc--; argv++;
  131.     if(strchr(argv[0], '?')) help();
  132. /*            open input file         */
  133.     if(argc && *argv[0]!='-')
  134.         {ifile=fopen(argv[0], "r");
  135.         if(ifile==0) {fprintf(stderr, "file %s not found\n", argv[0]); exit(1);}
  136.         strcpy(iname, argv[0]);
  137.         argc--; argv++;
  138.         }
  139. /*            open output file         */
  140.     if(argc && *argv[0]!='-')
  141.         {strcpy(oname, argv[0]);
  142.         argc--; argv++;
  143.         unlink(oname);
  144.         if((ofile=fopen(oname, "w"))==0)
  145.             {fprintf(stderr, "can\'t open output file %s", oname);
  146.             exit(1);
  147.             }        
  148.         }
  149.  
  150.     fprint_cmd(ofile, "; ifft %s\n");
  151.  
  152.     p_data[0]=-1;
  153.     i=1;        /* note x[0] and y[0] aren't defined */
  154.     while(i<ENTRIES)
  155.         {if(fgets(buf, BUFSIZE, ifile)==0) break;
  156.         t=buf;
  157.         while(*t && isspace(*t)) t++;
  158.         if(*t == 0) continue;        /* skip blank lines */
  159.         buf[strlen(buf)-1]=0;        /* zero out the line feed */
  160.         if(buf[0]==';')                /* skip comment */
  161.             {fprintf(ofile, "%s\n", buf); continue;
  162.             }
  163.         if(t = strchr(buf, ';')) *t=0;    /* zap same-line comment */
  164.         if(automatic_abscissas)
  165.             {xx=abscissa;
  166.             abscissa+=abscissa_step;
  167.             sscanf(buf, "%lf %lf", &yyr, &yyi);
  168.             }
  169.         else
  170.             {sscanf(buf, "%lf %lf %lf", &xx, &yyr, &yyi);
  171.             }
  172.         range_check(xx); range_check(yyr); range_check(yyi);
  173.         x[i]=xx; y[i+i-1]=yyr; y[i+i]=yyi; /* convert doubles to floats here */
  174.         s=buf;                      /* start looking for label */
  175.         nums=3;
  176.         if(automatic_abscissas) nums--;
  177.         while(nums--)                    /* skip the numbers */
  178.             {while(*s==' ')s++;
  179.             while(*s && (*s!=' '))s++;
  180.             }
  181.         while(*s==' ')s++;
  182.         if((length=strlen(s))&&(labels<MAXLABELS))
  183.             {if(*s=='\"')           /* label in quotes */
  184.                 {t=s+1;
  185.                 while(*t && (*t!='\"')) t++;
  186.                 t++;
  187.                 }
  188.             else                    /* label without quotes */
  189.                 {t=s;
  190.                 while(*t && (*t!=' '))t++;
  191.                 }
  192.             *t=0;
  193.             length=t-s;
  194.             p_data[labels]=i;
  195.             p_text[labels]=(char *)malloc(length+1);
  196.             if(p_text[labels]) strcpy(p_text[labels++], s);
  197.             }
  198.         i++;
  199.         }
  200.     n=i-1;
  201.     if(breaking && (!labels || p_data[labels-1]!=n-1))
  202.         {p_data[labels]=i-1;
  203.         if(p_text[labels]=(char *)malloc(1)) *p_text[labels++]=0;
  204.         }
  205. }
  206.  
  207. /* check whether number is too big for a float */
  208. range_check(x) double x;
  209. {    if(fabs(x)>BIGFLOAT)
  210.         {printf("input number too large: %f\n", x);
  211.         exit(1);
  212.         }
  213. }
  214.  
  215. int streq(a,b) char *a,*b;
  216. {    while(*a)
  217.         {if(*a!=*b)return 0;
  218.         a++; b++;
  219.         }
  220.     return (*b==0);
  221. }
  222.  
  223. /* get_parameter - process one command line option
  224.         (return # parameters used) */
  225. get_parameter(argc, argv) int argc; char **argv;
  226. {    int i;
  227.     if(streq(*argv, "-a"))
  228.         {i=get_double(argc, argv, 1, &abscissa_step, &abscissa, &abscissa);
  229.         abscissa_arguments=i-1;
  230.         automatic_abscissas=1;
  231.         return i;
  232.         }
  233.     else if(streq(*argv, "-z"))
  234.         {origin=0.;
  235.         i=get_double(argc, argv, 1, &origin, &fmax, &fmax);
  236.         return i;
  237.         }
  238.  
  239.     else gripe(argv);
  240. }
  241.  
  242. char *message[]={
  243. "ifft   version ", VERSION,
  244. " - calculate inverse fast Fourier transform \n",
  245. "                     of frequency domain function\n",
  246. "usage: ifft  [infile  [outfile]]  [options]\n",
  247. "options are:\n",
  248. "  -a  [step]     automatic abscissas \n",
  249. "  -z  val        add <val> to each abscissa value\n",
  250. /*
  251. "  -b             break input after each label, \n",
  252. "                 find separate transforms\n",
  253. "  -f  min [max]  minimum and maximum frequencies\n",
  254. "  -t  min [max]  minimum and maximum times\n",
  255. */
  256. 0
  257. };
  258.  
  259.  
  260. /*    name...
  261.         fft
  262.  
  263.     purpose...
  264.         perform forward or inverse FFT
  265.  
  266.     history...
  267.         May 84  translated from FORTRAN
  268.         12 Jun 84  diagnostic printouts shortened
  269. */
  270.  
  271.     int ip, n2, i, j, k, nu, nn;
  272.     int k1n2, l, nu1;
  273.     double arg, p, q, s, c;
  274.     double ninv, time, dt, freq, df;
  275.     double xplus, xminus, yplus, yminus;
  276.     double uplus, uminus, vplus, vminus;
  277.  
  278. /*    format frequency domain data so its
  279.     inverse fast Fourier transform can be calculated
  280.  
  281.     on input:
  282.     n is one more than a power of 2:  n == 3, 9, 17, 33, 65, ...
  283.     x[1...n]  has  f(0)  f(1)  f(2)  ...  f(n-1)
  284.     y[1...2n] has  Qr(0) Qi(0)   Qr(1) Qi(1)  ...  Qr(n-1) Qi(n-1)
  285.  
  286.     on output:
  287.     n has been increased:  n = 2*n-1
  288.     x[1...n] has  t(0)  t(1)  t(2)  ...  t(n-1)
  289.     y[1...n] has  q(0)  q(1)  q(2)  ...  q(n-1)
  290. */
  291. ifft(nnn, x, y) int *nnn; FLOAT x[], y[];
  292. {
  293.     n= *nnn;
  294.     if(n<= 0) {fprintf(stderr, "ifft: bad value of n: %d", n); return 1;}
  295.     df= x[2]-x[1]; dt= .5/x[n];
  296.     if((fabs(df*dt*(n-1)-.5)>.01) || (fabs(x[n]-x[n-1]-df)>.01*df))
  297.         {fprintf(stderr, "ifft: frequencies not equally spaced");
  298.         exit(1);
  299.         }
  300.     nn= 2; nu= 0;
  301.     while(++nu<= 20)
  302.         {if(nn==n-1)break;
  303.         nn= nn+nn;
  304.         }
  305.     if(nu>20)
  306.         {fprintf(stderr, "ifft: n isn\'t 1 more than a power of 2");
  307.         return n;
  308.         }
  309.     ip= 2; i= 0;
  310.     while(++i<= n) {x[i]= y[ip]; y[i]= y[ip-1]; ip=ip+2;}
  311.     n2= n/2+1; ip= n; arg= 0.;
  312.     n= n-1;   /* original n, minus one */
  313.     ninv= 3.14159265358979/n;
  314.     i= 0;
  315.     while(++i<= n2)
  316.         {/* printf("i= %d ", i); */
  317.         s= sin(arg); c= cos(arg);
  318.         yplus= y[i]+y[ip]; yminus= y[i]-y[ip];
  319.         xplus= x[i]+x[ip]; xminus= x[i]-x[ip];
  320.         p= s*yminus+c*xplus; q= s*xplus-c*yminus;
  321.         x[i]= q-xminus; x[ip]= q+xminus;
  322.         y[i]= yplus-p;  y[ip]= yplus+p;
  323.         --ip;
  324.         arg= arg+ninv;   /* faster than: arg=ninv*i; */
  325.         }
  326.  
  327.     fft2(y, x, n, nu);
  328.  
  329. /*    fprintf(stderr, "ifft: transform calculated\n");
  330.     i=0; while(++i<=10) {printf("%6d %15.6f %15.6f\n", i, y[i], x[i]);}
  331. */
  332.     ip=n;
  333.     while(ip)
  334.         {y[ip+ip]= -df*x[ip]; y[ip+ip-1]=df*y[ip];
  335.         --ip;
  336.         }
  337.     time=0.; i=0;
  338.     *nnn=n=n+n;  /* twice (original n, minus 1) */
  339.     while(++i<=n) {x[i]=time; time=time+dt;}
  340. }
  341. fft2(xreal, ximag, n, nu)
  342. FLOAT xreal[], ximag[];
  343. int n, nu;
  344. {    double treal, timag;
  345.     FLOAT *xr1, *xr2, *xi1, *xi2;
  346.  
  347.     n2=n/2; nu1=nu-1; k=0; l=0;
  348.     while(++l<=nu)
  349.         {while(k<n)
  350.             {p=ibitr(k>>nu1, nu);
  351.             arg=3.14159265358979*2*p/n;
  352.             c=cos(arg); s=sin(arg); i=0;
  353.             while(++i<=n2)
  354.                 {++k; k1n2=k+n2;
  355.                 /* initialize four pointers */
  356.                 xr1=xreal+k; xr2=xreal+k1n2;
  357.                 xi1=ximag+k; xi2=ximag+k1n2;
  358.                 treal= *xr2*c+ *xi2*s;
  359.                 timag= *xi2*c- *xr2*s;
  360.                 *xr2= *xr1-treal;
  361.                 *xi2= *xi1-timag;
  362.                 *xr1= *xr1+treal;
  363.                 *xi1= *xi1+timag;
  364.                 }
  365.             k=k+n2;
  366.             }
  367.         k=0; --nu1; n2=n2/2;
  368.         }
  369.     k=0;
  370.     while(++k<=n)
  371.         {i=ibitr(k-1, nu)+1;
  372.         if(i>k)
  373.             {treal=xreal[k];   timag=ximag[k];
  374.             xreal[k]=xreal[i]; ximag[k]=ximag[i];
  375.             xreal[i]=treal;    ximag[i]=timag;
  376.             }
  377.         }
  378. }
  379.  
  380. /* reverse the last nu bits of j */
  381. ibitr(j, nu) int j, nu;
  382. {    int ib;
  383.     ib=0;
  384.     while(nu--){ib=(ib<<1)+(j&1); j=j>>1;}
  385.     return ib;
  386. }
  387.  
  388.  
  389.